      SUBROUTINE MF_MDHNKL (z,h1,h2,h1p,h2p,expz,ngboth)

c***********************************************************************

c     This routine calculates the modified Hankel functions of order 1/3
c     and their derivatives.

c  Change History:
c     26 Oct 95     Changed to get the LOG unit from LWPC_LUN.CMN.

c*******************!***************************************************

      implicit complex (a-h,o-z)

c     LWPC parameters
      include      'lwpc_lun.cmn'

      complex       h1,h2,h1p,h2p,z,expz

      real          a(30),b(30),c(30),d(30),cap(30),zmag,ztlog
      complex       mpower,mterm

      real          part1(2),part2(2)
      equivalence  (part1,term4),(part2,sum4)

      data      a/9.3043671692922944819e-01,  3.1014557230974314911e+01,
     &            2.0676371487316209897e+02,  5.7434365242545027449e+02,
     &            8.7021765519007617234e+02,  8.2877871922864397320e+02,
     &            5.4168543740434246542e+02,  2.5794544638302022111e+02,
     &            9.3458495066311674231e+01,  2.6626351870744066662e+01,
     &            6.1210004300561072794e+00,  1.1592803844803233472e+00,
     &            1.8401275944132116616e-01,  2.4833030963741048003e-02,
     &            2.8842080097260218300e-03,  2.9133414239656786138e-04,
     &            2.5827494893312753646e-05,  2.0256858739853140063e-06,
     &            1.4155736366074870734e-07,  8.8695090013000443124e-09,
     &            5.0110220346327933889e-10,  2.5658074934115685526e-11,
     &            1.1961806496091228666e-12,  5.0988092481207283185e-14,
     &            1.9948392989517716388e-15,  7.1886100863126905797e-17,
     &            2.3938095525516785112e-18,  7.3883010881224645255e-20,
     &            2.1194208514407528762e-21,  5.6653858632471341093e-23/

      data      b/6.7829872514427588456e-01,  1.1304978752404598033e+01,
     &            5.3833232154307609704e+01,  1.1962940478735024376e+02,
     &            1.5337103177865415841e+02,  1.2780919314887846509e+02,
     &            7.4742218215718400631e+01,  3.2355938621523117060e+01,
     &            1.0785312873841039006e+01,  2.8532573740320209005e+00,
     &            6.1360373635097223595e-01,  1.0937678009821251966e-01,
     &            1.6422939954686564465e-02,  2.1055051223957133911e-03,
     &            2.3316778764072130571e-04,  2.2528288660939256561e-05,
     &            1.9156708045016374595e-06,  1.4446989475879618839e-07,
     &            9.7286124416697769730e-09,  5.8854279743918795891e-10,
     &            3.2160808603234314644e-11,  1.5952782045255116351e-12,
     &            7.2151886229105003778e-14,  2.9876557444763976717e-15,
     &            1.1368553061173507104e-16,  3.9889659863766691603e-18,
     &            1.2946984700995355913e-19,  3.8985199340546088228e-21,
     &            1.0920223904914870636e-22,  2.8527230681595795812e-24/

      data      c/4.6521835846461472410e-01,  6.2029114461948629822e+00,
     &            2.5845464359145262382e+01,  5.2213059311404570392e+01,
     &            6.2158403942148298012e+01,  4.8751689366390821897e+01,
     &            2.7084271870217123228e+01,  1.1215019407957400909e+01,
     &            3.5945575025504490022e+00,  9.1815006450841609147e-01,
     &            1.9128126343925335199e-01,  3.3122296699437809740e-02,
     &            4.8424410379295043444e-03,  6.0568368204246458321e-04,
     &            6.5550182039227768583e-05,  6.1985987743950608612e-06,
     &            5.1654989786625507119e-07,  3.8220488188402150986e-08,
     &            2.5278100653705126277e-09,  1.5033066103898380141e-10,
     &            8.0822936042464409157e-12,  3.9473961437101054471e-13,
     &            1.7590891906016512675e-14,  7.1814214762263778920e-16,
     &            2.6957287823672589641e-17,  9.3358572549515461865e-19,
     &            2.9922619406895981315e-20,  8.9015675760511620701e-22,
     &            2.4644428505125033375e-23,  6.3656020935361057409e-25/

      data      d/6.7829872514427588456e-01,  4.5219915009618392131e+01,
     &            3.7683262508015326776e+02,  1.1962940478735024344e+03,
     &            1.9938234131225040548e+03,  2.0449470903820554375e+03,
     &            1.4201021460986496090e+03,  7.1183064967350857463e+02,
     &            2.6963282184602597492e+02,  7.9891206472896585111e+01,
     &            1.9021715826880139294e+01,  3.7188105233392256682e+00,
     &            6.0764877832340288572e-01,  8.4220204895828535644e-02,
     &            1.0026214868551016149e-02,  1.0363012784032058021e-03,
     &            9.3867869420580235442e-05,  7.5124345274574017960e-06,
     &            5.3507368429183773360e-07,  3.4135482251472901638e-08,
     &            1.9618093247972931935e-09,  1.0209780508963274472e-10,
     &            4.8341763773500352579e-12,  2.0913590211334783723e-13,
     &            8.2990437346566602039e-15,  3.0316141496462685641e-16,
     &            1.0228117913786331176e-17,  3.1967863459247792364e-19,
     &            9.2821903191776400453e-21,  2.5103962999804300309e-22/

      data    cap/1.0416666666666666663e-01,  8.3550347222222222116e-02,
     &            1.2822657455632716019e-01,  2.9184902646414046315e-01,
     &            8.8162726744375764874e-01,  3.3214082818627675264e+00,
     &            1.4995762986862554546e+01,  7.8923013011586517530e+01,
     &            4.7445153886826431887e+02,  3.2074900908906619004e+03,
     &            2.4086549640874004605e+04,  1.9892311916950979121e+05,
     &            1.7919020077753438063e+06,  1.7484377180034121023e+07,
     &            1.8370737967633072978e+08,  2.0679040329451551508e+09,
     &            2.4827519375935888472e+10,  3.1669454981734887315e+11,
     &            4.2771126865134715582e+12,  6.0971132411392560749e+13,
     &            9.1486942234356396792e+14,  1.4413525170009350101e+16,
     &            2.3788844395175757942e+17,  4.1046081600946921885e+18,
     &            7.3900049415704853993e+19,  1.3859220004603943141e+21,
     &            2.7030825930275761623e+22,  5.4747478619645573335e+23,
     &            1.1498937014386333524e+25,  2.5014180692753603969e+26/

      data        root3/(1.73205080756888,0.)/
     &            alpha/(8.53667218838951e-1,0.)/


      zz=z
      zmag=ABS(zz)
      zpower=(1.,0.)
      if (zmag .le. 4.2) then                 ! Power series expansion

         if (zmag .le. 1.e-10) then

            sum1=a(1)
            sum2=b(1)
            sum3=c(1)
            sum4=d(1)
         else

            zterm=-zz**3/(200.,0.)
            ztlog=LOG(zterm)
            sum1=0.
            sum2=0.
            sum3=0.
            sum4=0.
            m=1
            last=0
            do while (m .le. 30 .and. last .eq. 0)

               sum1 =sum1+a(m)*zpower
               sum2 =sum2+b(m)*zpower
               sum3 =sum3+c(m)*zpower
               sum4 =sum4+d(m)*zpower
               term4=     d(m)*zpower
               if (ABS(part1(1)) .le. 1.e-5*ABS(part2(1)) .and.
     &             ABS(part1(2)) .le. 1.e-5*ABS(part2(2))) then

                  last=1
               else
     &         if (REAL(LOG(zpower))+ztlog .lt. -30.) then

                  last=1
               else

                  zpower=zpower*zterm
               end if
               m=m+1
            end do
         end if

Change: 17 Oct 93; The following 2 statements were too complicated for WATCOM's
c                  compiler, version 9.5.

c         gm2f=(0.,1.)*(zz*sum2-(2.,0.)*sum1)/root3
c         gpmfp=(0.,1.)*(sum4+(2.,0.)*zz*zz*sum3)/root3

         gm2f=zz*sum2-(2.,0.)*sum1
         gm2f=(0.,1.)*gm2f/root3

         gpmfp=sum4+(2.,0.)*zz*zz*sum3
         gpmfp=(0.,1.)*gpmfp/root3

Change: 17 Oct 93; End

         zh1=zz*sum2+gm2f
         zh2=zh1-(2.,0.)*gm2f
         zh1p=sum4+gpmfp
         zh2p=zh1p-(2.,0.)*gpmfp
         zexp=(0.,0.)
      else                                        ! Asymptotic expansion

         mpower=(1.,0.)
         sum1=(1.,0.)
         sum2=(1.,0.)
         sum3=(0.,0.)
         sum4=(0.,0.)
         rootz=SQRT(zz)
         rootz_cubed=rootz*zz
         zterm=(0.,1.)/rootz_cubed
         mterm=-zterm
         dm=(0.,0.)
         term3=(1.,0.)
         last=0
         m=1
         do while (last .eq. 0 .and. m .le. 30)

            zpower=zpower*zterm
            mpower=mpower*mterm
            dm=dm+(1.,0.)
            term1=cap(m)*zpower
            term2=cap(m)*mpower
            if (ABS(term2/term3) .ge. 1.) last=1
            sum1 =sum1+term1
            sum2 =sum2+term2
            sum3 =sum3+term1*dm
            sum4 =sum4+term2*dm
            term4=     term2*dm
            if (ABS(part1(1)) .le. 1.e-5*ABS(part2(1)) .and.
     &          ABS(part1(2)) .le. 1.e-5*ABS(part2(2))) last=1
            term3=term2
            m=m+1
         end do
         sum3=sum3*zterm*(-1.5,0.)/zz
         sum4=sum4*zterm*(-1.5,0.)/zz

         term1=((-0.25,0.)-(0.,1.)*rootz_cubed)/zz
         term2=((-0.25,0.)+(0.,1.)*rootz_cubed)/zz

         zh1 =sum2
         zh1p=sum2*term2+sum4
         zh2 =sum1
         zh2p=sum1*term1+sum3

c        zexp=-i*2/3*z**(3/2)+i*pi*5/12
         zexp=(0.,-0.666666666666667)*rootz_cubed
     &       +(0., 1.308996938995747)

         if ( REAL(zz) .lt. 0.) then

            exp1=EXP(zexp)
c           exp2=EXP(zexp-i*pi*4/3)
            exp2=EXP(zexp+(0.,-4.188790207))

            if (AIMAG(zz) .ge. 0.) then

               zh2 =zh2 *exp1
               zh2p=zh2p*exp1

               if (ngboth .eq. 0) then

                  zh2 =zh2 +zh1 /exp2
                  zh2p=zh2p+zh1p/exp2
               end if

               zh1 =zh1 /exp1
               zh1p=zh1p/exp1
            else

               th2 =-zh1 /exp2
               th2p=-zh1p/exp2

               zh1 =zh1 /exp1+zh2 *exp2
               zh1p=zh1p/exp1+zh2p*exp2

               if (ngboth .eq. 0) then

                  zh2 =zh2 *exp1
                  zh2p=zh2p*exp1
               else

                  zh2 =th2
                  zh2p=th2p
               end if
            end if
            zexp=(0.,0.)
         end if

         zterm=alpha/SQRT(rootz)
         zh1  =zterm*zh1
         zh2  =zterm*zh2
         zh1p =zterm*zh1p
         zh2p =zterm*zh2p
      end if
cxxcDEBUG
cxxc     Calculate Wronskian as partial check on validity;
cxxc     this test is not valid when h1, h2, h1', h2' are of the same
cxxc     order of magnitude.
cxx      term4=zh1*zh2p-zh1p*zh2
cxx      if (ABS(part1(1))                   .gt. 1.e-4 .or.
cxx     &    ABS(part1(2)+1.457495441040461) .gt. 1.e-4)
cxx     &   write(lwpcLOG_lun,
cxx     &       '(''WARNING: [MF_MDHNKL]:''/(a,1p2e15.6))')
cxx     &         'z',z,'W',sum4

c     Return values in single precision
      h1  =zh1
      h2  =zh2
      h1p =zh1p
      h2p =zh2p
      expz=zexp
      RETURN
      END      ! MF_MDHNKL
